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We apply the self-consistent harmonic approximation (SCHA) to study static and dynamic 
properties of the two-dimensional classical Heisenberg model with easy-axis anisotropy. The 
static properties obtained are magnetization and spin wave energy as functions of temper- 
ature, and the critical temperature as a function of the easy-axis anisotropy. We also 
calculate the dynamic correlation functions using the SCHA renormalized spin wave energy. 
Our analytical results, for both static properties and dynamic correlation functions, are 
compared to numerical simulation data combining cluster-Monte Carlo algorithms and Spin 
Dynamics. The comparison allows us to conclude that far below the transition temperature, 
where the SCHA is valid, spin waves are responsible for all relevant features observed in the 
numerical simulation data; topological excitations do not seem to contribute appreciably. 
For temperatures closer to the transition temperature, there are differences between the 
dynamic correlation functions from SCHA theory and Spin Dynamics; these may be due to 
the presence of domain walls and solitons. 

PACS numbers: 75.10.Hk; 75.40.Cx 



I. INTRODUCTION 

Low-dimensional magnets have been extensively inves- 
tigated by many theorists and experimentalists in the 
last three decades. More recently, the interest on the 
properties of two-dimensional(2D) Heisenberg magnets 
has been greatly revived since the discovery of high-Tc 
superconductivity: it is now well knownEJ that the un- 
doped, insulating La2Cu04 has a quasi-two-dimensional 
antiferromagnetic behavior. However, most quasi-two- 
dimensional magnetic real materials exhibit some kind 
of anisotropy: the anisotropic properties often arise not 
so much from an anisotropy in the interaction mechanism 
(which can be wholly isotropic) but from other sources, 
such as the presence of a crystal field that couples the 
spins to a certain direction in the crystal. Then, at least 
from a theoretical point of view, a large amount of mag- 
netic materials fits (under certain circumstances like tem- 
perature range) into one of the two groups: easy-plane 
or easy-axis models. Easy-plane 2D magnets have de- 
served a lot of attention due to their possibility of show* 
ing the topological Kosterlitz-Thouless phase transition.u 
The interest devoted to easy-axis magnetic systems has 
been considerably smaller, specially concerning the study 
of its dynamical properties. It is our aim to address to 
this topic in this paper. 

It must be emphasized that, although we shall be con- 
cerned only with magnetic systems in this paper, many 
of the magnetic Hamiltonians also allow for an interpre- 
tation other than a magnetic one. Most physical prob- 
lems concerning mutually interacting elements that form 
a spatial array can be mapped into a magnetic Hamil- 



tonian by describing it within a pseudo spin formalism. 
The advantage of studying a general physical problem 
in its magnetic form is clearly that in magnetism sev- 
eral experimental techniques are amilable to study the 
fundamental properties of a system.S 

The analysis of the general Ising-Heisenberg model is of 
interest because, from the experimental point of view, the 
presence of some degree of anisotropy in the interaction 
is to be expected in nearly all cases. In addition, recently 
there has been a growing interest in the study of topolog- 
ical excitations in the classical two-dimensional easy-axis 
modcl.Q Having finite excitation energy, the population of 
topological objects should be quite small at low tempera- 
tures. Therefore, before taking into account the effect of 
topological excitations (solitons or similar objects) on the 
thermodynamics and dynamics of a system, we should 
consider the contribution of anisotropic spin waves. So 
we might ask: can spin waves explain experimental data 
or, in the absence of experiments, computer simulation 
data at low temperatures? This is the spirit and aim of 
this paper. 

Here we consider the classical Heisenberg ferromag- 
net in two dimensions (2D) with easy-axis exchange 
anisotropy 



H — —J 2_^ S„ • Sn+a — K 2_^ S^S^+a 
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where the summations run over all distinct couples of spin 
sites n and its nearest neighbors a. As the anisotropy pa- 
rameter K ranges from to cx), we go from the isotropic 
Heisenberg model to an Ising like model in which the 
spins tend to be confined along the ±z— direction. How- 



ever, the resemblance to the Ising behavior can only hold 
for T < AT: we find that, for Hamiltonian ^, Tc ^ K 
for large K. This contrasts with T^ w 2.27 K for the 2D 
single-component Ising model. 

In addition to the usual domain walls we expect that 
there can be localized soliton-like excitations that can 
connect a small circular domain of positive S^ to a sur- 
rounding region with negative 5^. A spatial "width" 
of these objects (bubbles or droplets) can be estimated 
as approximately \J J jK . For intermediate values of 
yJj/K, i.e., between a lattice constant and the system 
size, these excitations can be important on a finite dis- 
crete system. These objects can also have a topological 
charge or winding number of the spin field. There was 
some indication in earlier Monte Carlo (MC) simulationsu 
that they may play a role in the phase transition in 
this model; their density was found to increase strongly 
passing through the transition temperature. However, 
in a continuum static description they are found to be 
energetically unstable, according to the Derrick-Hobart 
theoremjj Thus it makes sense to investigate whether it is 
necessary to be aware of their presence in static and dy- 
namic properties of this model, or whether a description 
based on anisotropic spin waves is sufficient. 

To this end we study the low temperature thermody- 
namics and dynamics of this model using a self consistent 
harmonic approximation theory (SCHA) to treat spin 
waves. As is well known, the SCHA is a reasonable ap- 
proximation to calculate the transition temperature and 
low temperature (T < Tc) properties of a system but 
it is of limited value in estimating critical properties. 
Therefore, in our work, we did not attempt to do any 
calculation for critical exponents and related aspects of a 
phase transition. We compare the predictions of SCHA 
theory to numerical simulations on several L x L square 
lattices (L = 16, 32, 64, 128) using Monte-Carlo and spin- 
dynamics (SD) simulations, which include effects due to 
all thermodynamically allowed excitations. We present 
the thermodynamic results in Section O, and in Section 
III, the calculation of the dynamical correlation function. 
The simulation procedures are discussed in Sections II B 



and [II B, and their comparison with the SCHA theo- 
retical calculations is given in Section [V. Finally, our 



conclusions are given in Section ^ 



spin waves, being characterized by simple temperature- 
dependent renormalization factors for the unperturbed 
spin wave energy. 

We start by writing the spin components using the 
Dyson-Maleev representation of spin operators 
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5ji — O flji^n 



where ajj and On are the Bose spin operators on site n. 
The harmonic spin wave Hamiltonian obtained from (|l|) 
is given by 



iJo = ^ Wqajjflq 



(3) 



where a], and Oq are the Fourier transforms of aj^ and 
respectively, and 



^q = 4J5[1 - 7(q)] + 4A'5 



(4) 



with 7(q) = ^[cosqx + cosqj^]. The spin wave approxi- 
mation will be reasonable when (aJjOn) <C 5, so it ought 
to be fairly good for anisotropics satisfying the relation 
T < 4:KS^. 

Now we simplify the general model by reducing 
Hamiltonian (|l|) to an effective harmonic problem with 
the effect of anharmonicity embodied in temperature- 
dependent renormalized parameters. This means that 
the couplings of the model are replaced by quadratic in- 
teractions whose strength is then optimized. ^Details of 
this method may be found in the literatureQu and we 
give here only an outline of those steps pertinent to our 
present calculation. 

We assume as effective Hamiltonian the appropriate 
form for a noninteracting gas of Bose excitations 
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The spin wave energy is obtained by a variational proce- 
dure based on the inequality for the Free energy F 



II. STATIC PROPERTIES 

A. Self-Consistent Harmonic Approximation 

Since its original derivation by Bloch,Q the self con- 
sistent harmonic approximation has been found to ac- 
count for the low temperature dependence of various 
properties of several magnetic insulators, which sgep. 
to be fairly well-described by the Heisenberg model.ETtS 
Its usefulness stems mainly from the way it takes into 
account a substantial part of the interactions among 
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where the brackets indicate the thermal average. Traces 
should be taken only over the physical states, that is, 
states with no more than 25 spin deviations on a single 
site. The minimization of (^ with respect to Eq deter- 
mines the spin wave energies.- We obtain, in the classical 
limit, following Rastelli et al,Ll 



i?q(r) = 4J5(l-7(q)) l-P{T) + f,{T) 
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where TV is the number of sites. Eqs. (R), (ph, and (H) 
are coupled equations which we solved self consistently 
by an iterative method. These coupled equations have 
a double-valued solution below Tc and no real solution 
above Tc- this is the typical behavior for self consistent 
harmonic approximations and allows for easy determina- 
tion of Tc- The lower branch (for T < Tc) has an unphys- 
ical temperature dependence and may be discarded as a 
spurious mathematical solution that is physically unsta- 
ble. In Figure Q^ the spin wave energy for K/J = 0.05 
is given for two temperatures well below T^ ~ 0.75J: 
T = 0.3J, and T = 0.6J. The circles and stars shown in 
Figure [l| were taken from our numerical simulation data 
(to be described in Section IIIB). As can be seen, the 
comparison between the SCHA and numerical results is 
remarkably good: the SCHA describes well the decrease 
of energy with increasing temperature (and, also, the en- 
ergy dependence with the wavevector). 

The reduced spontaneous magnetization along the z- 
axis is given by 
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In Figure g, we present results obtained from Eq. ( [lO| ) for 
K/J = 0.05 and compare to our Monte Carlo (MC) data 
(obtained as described in Section II B ). The slight over- 
estimate of Tc from SCHA clearly is due to the fact that 
it does not include all possible modes of fluctuations, that 
are included in the MC calculations. The SCHA theory 
has no built-in requirement to make the magnetization 
null at the critical temperature Tc^ and consequently, we 
find, as Figure || shows, a nonzero value for Mz(Tc). This 
is typical of SCHA approaches: the theory applies only 
below Tc, and for temperatures T > Tc the magnetization 
is taken as zero, implying a discontinuous jump at Tc- In 
fact, the scaling of the MC data for Mz{T) with system 
size L strongly suggests the presence of a discontinuous 
jump. 



B. Monte Carlo 

In order to evaluate the accuracy of the above the- 
ory, we calculated Tc and the magnetization and other 
thermodynamic quantities using a hybrid classical Monte 
Carlo approach on periodic Lx L square lattices. We ap- 
plied a combination of Metropolis single-spin moves and 



over-relaxation moves that modify all three spin compp=. 
nents, and in addition, Wolff single-cluster operationsEJ 
that modify only the S^ components. The over- 
relaxation and cluster moves are necessary to avoid crit- 
ical slowing down near Tc, which is tending to freeze 
the S^ components. The single spin and over-relaxation 
moves are standard, here we give only a few details about 
the cluster algorithm. In the Wolff single-cluster algo- 
rithm, the cluster-flip operation we used only reverses 
the sign of S^ for all sites that have been included into 
the clustcp^ This is reminiscent of the Swendsen-Wang 
algorithmlla for Ising models, but we only build one clus- 
ter at a time as in the Wolff algorithm. The cluster moves 
cannot be used alone because they do not change the 
magnitudes of S^ spin components. 

A cluster is built up starting from a randomly cho- 
sen seed site n, immediately inverting its 5*^ component: 
S^ — > — S'n, and then including neighboring sites n + a 
with a probability. 
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Here A_En,n-i-a is the energy change involved if site n -|- a 
is not flipped: 



A£;„.„+a = 2(J + i^)5^5^ 



(12) 



Note that in this formula site n was already included 
into the cluster and S^ was already inverted. Eq. ( |ll| ) 
represents the cluster growth as essentially a sequence of 
Metropolis decisions, according to whether A_En^„+a is 
less than or greater than zero. Newly included sites then 
have their neighbors tested for inclusion until the cluster 
is done growing, at which point all included sites have 
already been modified. 

We define one cluster sweep as building enough sin- 
gle clusters until the number of sites included into clus- 
ters is one quarter of the total number of sites in the 
system. Then we defined one hybrid Monte Carlo step 
as one over-relaxation sweep followed by one Metropo- 
lis single spin sweep followed by one Wolff cluster sweep. 
Equilibrium data shown here are averages over 10^ to 
4 X 10^ Monte Carlo steps. The critical temperature 
was determined from the change in the distribution of 
z-component of total magnetization, which is-easily char- 
acterized by Binder's fourth cumulant ratio,t3 
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The crossing point of curves of Ul (T) for different system 
sizes gives a good estimate of Tc- All calculations were 
made for square lattices of size L x L, using unit spins 
5 = 1 and fixing J — 1 while allowing K to be varied. 



C. Static Results 

The critical temperature from the SCHA as a function 
of anisotropy parameter K/J is shown in Figure ^ and 



compared with numerical MC estimates for a set of spe- 
cific values of K/ J ranging from 0.05 to 10.0. Generally, 
the SCHA overestimates T^ when compared to the MC 
results, because it does not fully take into account all pos- 
sible fluctuations that are taking part in the transition. 
Notice that, as K increases, the dependence of Tc on K/ J 
becomes linear. For K/ J 3> 1, we recover a continuous 
spin Ising Haniiltonian: Eq. (yj) can be approximated as 
i/ w J(l + K/J)S^S^+^ = KS^S^+^. Figure | shows 
that, for K/J > 1.0, the results follow a straight line with 
slope « 1.0. We remark again that, strictly speaking, the 
analogy between Hamiltonian (|l|) and the continuous spin 
Ising model can only be expected to hold for T <^ K. For 
moderate and high temperatures, model (|l|) still exhibits 
the full entropy effects of a three-component spin model, 
resulting in a much lower transition temperature than a 
one-component Ising model. For this reason we cannot 
expect to compare our results to the ones obtained for 
the usual 2D Ising model. 

Some of the drawbacks of the SCHA are well known: 
(i) it does not take into account strong coupling effects 
which are important at high temperatures and at short 
wavelengths; (ii) it also neglects the kinematical interac- 
tion and gives a first-order phase transition to the para- 
magnetic phase (where the true phase transition should 
be of second order). Notwithstanding this, we see that 
the theory gives results which compare quite well with 
the MC data we obtained. 

This good agreement cannot be used to conclude that 
solitons do not have an important contribution to the 
properties of our model. As is well known, in the one 
dimensional easy- axis ferromagnet, the soliton connects 
two distinct ground-states and has, therefore, a global ef- 
fect in the system.Eil As a consequence, a pure spin wave 
calculation does not predict correctly all thermodynamic 
quantities. For instance, spin waves give a linear be- 
havior with temperature T (for T — > 0) for the corre- 
lation length, while the soliton model predicts correctly 
an exponential behavior. In two dimensions, however, 
the soliton has only a local effect and its contribution to 
thermodynamic quantities should be small. The reason- 
able agreement of the SCHA calculation with the Monte 
Carlo data is not, per se, an indication that we should 
rule out topological effects. The signature of topological 
solitons is best analyzed in the dynamics where it should 
manifest as a central peak. This topic will be discussed 
in the following Sections. 
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The averages are readily evaluated, and give the xx- and 
yy- dynamical correlation functions: 



(^q(i)^-q> = (^^(i)^^q> 



(16) 
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where Uq is the Bose occupation number and E^q is the 
self-consistent spinwave frequency of Eq. (|^). Eq. ( |lm 
leads to pure spin wave peaks for the spectral functioiilj 
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The comparison between these SCHA results for S' ^^ and 
those obtained by spin dynamics simulation (Sec. IIIB) 
are shown in Figures n^ and H. These are discussed in 
more detail in Sec. |^.B below. 

The zz- correlation contains, in addition to the Bragg 
scattering at q = 0, a term describing correlations in the 
fluctuations of the spin's z component. Evaluating the 
averages, we obtain 
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where 



^ - ^a+k - -Ea-i, 



(18) 



(19) 



Eq. (|l^) corresponds to the various two-spin-wave scat- 
tering terms. It is interesting to notice that, to this order, 
only difference processes contribute to the dynamics. The 
time Fourier transform of ( |l8[) , together with the theripo- 
dynamic limit L ^ oo, gives us the response functioElld 



III. DYNAMIC CORRELATION FUNCTIONS 



A. SCHA 
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(20) 



From Hamiltonian M) and using Eqs. (0), we obtain 
the time dependent correlation functions 

{S^it)S\) = f (K(i)+<(i)][«-q + «U]> , 



Using the delta functions we obtain integrals on the con- 
tours C* defined by w = ±0: the first integral is 
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(21) 



where dli^ is the contour element and | Vfef^ | designates 
the Jacobian of the involved transformation. There is 
a singularity in these integrands for every minimum or 
maximum of Q and the spectrum is in general quite com- 
plicated. We emphasize that, here, we have also used the 
self-consistent result obtained in Section |l| for the spin 
wave energies. Results for the spectral functions obtained 
from (pol ) will be compared with those from MC-SD sim- 
ulations below in Section IV. 



B. Spin Dynamics Simulations 

The spin dynamics simulation is standardJl3il3 Here 
we summarize the method and describe the particular nu- 
merical parameters used. For a given temperature, a set 
of 200 initial states was taken from the Monte Carlo simu- 
lation to serve as initial conditions for the spin-dynamics 
time integration. The nonlinear equations of motion as- 
sociated with Hamiltonian p) are 



dt 
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(22) 



hand, for wavevectors q = ^(2m -I- 1,0), there is an in- 
tensity minimum in S'^^(q, cj) at w — > 0. In order to see if 
the SCHA theory could explain this interesting result we 
re-started our calculation from ([l8|), restricting the sums 
to the discrete wavevectors k = =^{m, n) of each lattice. 
Also, the time integration for 5(q, w) was performed for 
a finite time interval tmax , 
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where tmax was taken to be the same as the integration 
time (430/ J) used in our simulations. Eq. ( [18|) is mod- 
ified for a finite time interval, and a complete analysis 
leads to 
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where J is the diagonal matrix of exchange couplings, 

(23) 
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J= \ J 
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These were integrated forward in time using a standard 
fourth order Runge-Kutta scheme with time step h = 
0.035/ J (for small K/J). By saving data for time Fourier 
transforms at intervals dt = 6h, allows for measuring 
S{q,uj) out to Wmax — ^TT / dt « 30 J. We saved a total 
of Nt = 2^^ samples in time, integrating out to final 
time iinax = Ntdt « 430/ J, giving a frequency resolution 
of 6llj — 27r/tniax ~ 0.015 J. The space and time Fourier- 
transformed spin-spin correlations were averaged over the 
200 initial states to get S{q,uj), for both in-plane and 
out-of-plane spin components. 



IV. DYNAMIC CORRELATIONS: RESULTS 

A. Small Lattices (L < 64) 

At low temperatures T ^ Tc, we especially expect that 
the SCHA should give good agreement with the spin- 
dynamics simulation. We had noticed, however, that 
spin dynamics for small lattices gives an interesting set 
of unevenly spaced peaks in S'^^(q, lj), in contrast to one 
sharp peak at the spinwave frequency in S^^ (q, uj) , and 
also in contrast to the smooth behavior predicted for 
S'^^(q, w) by Eq. ([20|). Also, an intensity maximum in 
S^^ (q, uj) for w ^ is present for wavevectors of the 
form q — -^(2?Ti, 0), where m is an integer. On the other 



The expression can be thought to represent S{q,uj) as a 
sum over a set of narrow peaks of width approximately 
'2/tmaxi centered at frequencies il, determined by choos- 
ing k such that both § + k and ^ — k in Eq. ( [l9|) are 
allowed discrete wavevectors. Besides restricting the sum 
in (p5| ) to the discrete set of lattice wavevectors, the fi- 
nite time integration tmax implies discrete frequency in- 
crements Suj = Itx jtmax ~ 0.015J, the same as in our 
spin dynamics simulation. 

Examination of (p^ ) and ( psj ) allows us to conclude 
that a nonzero intensity in S^^ (q, w — > 0) can exist 
for all wavevectors and not only for those of the form 
q — ^(2771,0). However, a little consideration shows 
that if q/2 does not fall on a reciprocal lattice vector, 
then it is impossible to choose a value of k in Eq. (jlj) to 
give 51 = 0. Therefore, for wavevectors q = ^{2m+l, 0), 
none of the multiple peaks in (p5| ) will be centered at 
zero frequency, and 5'^^(q, w — *■ 0) is a local minimum. 
Although no peak is centered at zero, the tails can con- 
tribute there. On the other hand, for wavevectors such 
as q = ^(2m, 0), and q = ^{iri, m), we see that ^ falls 
on a highly symmetric point in the reciprocal lattice, and 
it is always possible to choose k to get J7 = in Eq. (|lS|). 
Then for these cases, there is a peak at zero frequency, 
and S^^{q,uj ^ 0) is a local maximum. 

The overall behavior of S^^{q,uj) with the lattice size 
obtained either by numerical simulation (Fig. o) or by the 
calculation of Eq. (|2^) (Fig. y) agree very well. In or- 
der to make this comparison, because the spin-dynamics 
simulations are purely classical, it is necessary to replace 
all factors of (1 -I- Uq) in the SCHA expressions by n^. 
Also, these occupation numbers were evaluated by their 
classical limit, riq = T/Eq, consistently with the static 



calculations in Sec. y. Figures |^ and ^ were obtained 
for K/J = 0.05 (where T^ « 0.75 J), q = (0.393,0), and 
T = 0.3J for lattice sizes L = 16, 32, and 64. 

Comparing the several peaks shown in Figures |5| and y, 
for a specific value of L, we see that they are positioned 
around the same frequencies. The important feature is 
that as the system size is increased, the spacing between 
the multiple peaks in S^^{q,uj) becomes smaller as -^. In 
addition, for a longer time integration tmax, the widths 
of the peaks will be narrower, and therefore they will be- 
come more distinct. As far as we aware, this strong finite- 
size effect in low-temperature spin-dynamics simulations 
is a feature that has been previously ignored. It is very 
likely, however, that it appears in any related models. 
For example, finite-size effects most likely explain simi- 
lar low-temperature multiple-peak features that have ap- 
peared in S{q, oj) calculatecLJ[p,the 2D Heisenberg model 
with easy-plane anisotropy.tSEj 



B. Large Lattices (L > 128) 

The SCHA calculation [Eq. (0)] and the MC-SD 
simulations both give single narrow spinwave peaks in 
S^^ (q, w) , regardless of lattice size. The MC-SD peak po- 
sitions for L = 128 have been compared with the SHCA 
results in Fig. nl and agree very well for the temperatures 
studied. The SCHA theory gives peaks of zero width, 
thus it makes sense to compare the integrated intensities 
for the positive frequency peak, I^^ = J^ °° dwS^^{q, lu). 
These are shown in Fig. |j, where the MC-SD results 
are compared to those obtained from Eq. O, J^^ = 
S'nq/(87r^). For the lower temperature, T = 0.3 J, there 
is very good agreement. The good low-T agreement, with 
no adjusted parameters, shows that the approximations 
made in the SHCA theory are reasonable where we ex- 
pect this simple theory to work. For T — 0.6J, how- 
ever, the MC-SD result is suppressed compared to the 
SCHA prediction. Currently we cannot say whether this 
suppression should be better described by spinwave in- 
teraction terms or possibly by nonlinear excitations such 
as solitons or domain walls. Clearly, both effects could 
become more important as the critical temperature is ap- 
proached. 

For S'^^(q, uj), the widths of the multiple peaks are de- 
termined both by the intrinsic width due to temperature, 
and the width 2T:/tmax inherent in the spin-dynamics 
simulation. For larger lattices, or higher temperatures, 
the spacing of the multiple peaks in S'^^(q, w) becomes 
smaller than their measured widths, the peaks merge and 
the curve is much smoother. Thus in our simulations the 
finite-size effects are quite well smoothed out for lattices 
L > 128 and/or for high temperatures. In Figures |7-10, 
for K = 0.05 J, several q values, and L = 128 we see that 
the simulation data for the higher temperature T — 0.6 J 
are smooth while the data for T = 0.3 J still show sharp 
peaks. 



Using the "discrete" equation (|2^) for obtaining 
5^^(q, w) for lattice size L — 128 we do not get rid of 
the multiple peak structure even for T = 0.6J. This 
can be seen in Fig. |l^ where the three types of calcula- 
tions — numerical simulation, discrete summation (Eq) 
and continuum limit (20) — we used to obtain ^^^(q, oj) 
are shown for K = 0.05, T = 0.6 J and q = (1.03,0). 
Typically, the discrete SCHA summation results in an 
S^^ (q, oj) curve with very strong multiple peak struc- 
ture. In order to smooth out the structure obtained 
from (p^) it is necessary to consider much larger lattices 
{L > 500). It is natural to expect that it is more diffi- 
cult to smooth out the spectra obtained by ( p5| ) than the 
one obtained via spin dynamic simulation. Clearly the 
MC-SD calculation contains more fluctuations and there- 
fore greater peak widths, especially as T approaches Tc, 
whereas in expression ( |25|) all spinwave peaks have very 
narrow widths determined only by the integration time. 
Instead of trying to smooth the SCHA spectra by consid- 
ering larger and larger lattices for the calculation of (p3) 
— which requires extra computational effort — we can 
go to the continuum approximation limit built in (EQ). In 
fact, most real systems contain a large number N of spins 
{N ^ oo) and effects due to the discreteness of the lat- 
tice are not important. These macroscopic systems will 
be better represented by the continuous approximation 
built in (|^). Figures ^|l~ 
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show the spectral functions 
obtained by numerically evaluating (|2y) for K = 0.05, 
T — 0.3J, and T = 0.6J, for the following wavevec- 
tors: q= (0.393,0), (1.473,0), (2.50,0) and (1.03,1.03). 
These are compared with the corresponding MC-SD cal- 
culations for 128 X 128 lattices. 

Obviously, considering the dynamical simulation, it is 
not possible to go to the N —^ oo limit: the computa- 
tional cost in simulations increases tremendously with N. 
Nevertheless, we can remark on interesting features con- 
cerning the results obtained from the SCHA calculation 
and from numerical simulation procedures. First, the fre- 
quency width of the region in which S'^^(q, uj) has appre- 
ciable intensity does not depend on the lattice size and 
on the kind of calculation performed to obtain S^^ (q, w) . 
This can be observed in Figs. ^, ^ and ^ which corre- 
spond to the three different ways we have used to obtain 
the spectral function for different lattice sizes but for the 
same wavevector q = (0.393, 0). 

In Fig. Il^ we show the comparison of the width Aw of 
the obtained spectral functions in the whole | q | range 
for wavevectors like q = {q, 0): the data were obtained for 
K = 0.05 J and T — 0.3 J. The comparison is remarkably 
good (a similar agreement is obtained for T — 0.6J). 
We see that, for small | q |, the width Alu increases 
linearly with the wavevector. A trivial analysis of (|2y) 
leads us to the conclusion that Aw must be related to 
the maximum value fi can have for each q. From ( |l9| ) 
we easily obtain that i^lmax = B{T) sin \ q \ /2 where 
B{T) = €JS[1 - (i{T) + r]{T){l + K)] and e = 1 for 
q = (9,0) wavevectors and e = 2 for q = {qx,qx)- For 
comparison, we show, in Fig. 03 a curve (dashed line) 



corresponding to flmax- 

Second, the SCHA curves corresponding to wavevec- 
tors of the form q = {q, 0) or (0, q) ( Figures ^^ ) show a 
sharp peak at higher frequencies, just before the spectral 
function vanishes. For wavevectors hke q = {q,q), this 
sharp peak is only observed near the point (tt, tt), and not 
for smaller |q|, as can be seen in Fig. 1^. The appearance 
or not of these peaks in the SCHA calculation depends 
on the behavior of the density of states | Vfi |^^ in Eq. 
(pi|). Figures |l^ and |lj show contours of 0(q, k) in the 
k-plane, for q ~ (0,0.393) and q = (1.03,1.03), respec- 
tively, for T = 0.3 J . For q in the (10) or (01) directions, 
the contours are straight lines (Fig. Ol). They become 
very widely spaced near k^ — 7r/2, or near k^ — 7i'/2, 
where ft approaches ftmax, and | VJ7 |~^ becomes very 
large along the whole straight contour. The integration 
along the contour in Eq. (|2l| ) then leads to the sharp 
peaks at w ^ ^max seen in ^^^(q, w). For q along the 
(11) direction, the contours are curves (Fig. |lj). For 
moderate values of | q |, the higher contours (near Qmax) 
approximate small circles, having limited total length and 
thus creating no sharp peak in 5'^^(q, lu). Only for q very 
close to the point (tt, tt) is the effect due to the divergence 
of I Vri |~^ more important than the contour length, and 
there a sharp peak at tj — > flmax does occur. 

It is interesting to notice that this peak is also present 
in the data obtained from the discrete spin wave calcu- 
lation (Figure |ll|) although the density of states does 
not appear explicitly in (p5f). Nevertheless, the two spin 
wave calculations must, in fact, give similar results be- 
cause (EQ) or ( Blf) c orrespond to the tmax —* oo, and 
iV ^ oo limits of (|25|). For small wavevectors, this sharp 
high frequency peak is not seen in the simulation data 
suggesting that inclusion of higher order terms in the 
spin wave theory would probably lead to the attenuation 
of this peak in the SCHA results. As the wavevector 
I q I increases, a lateral shoulder develops in the spectra 
obtained from both numerical simulation and SCHA cal- 
culations; it is already well defined for q ^ 0.50. For very 
large wavevectors, as in Figure H the lateral shoulder for 
the MC-SD data occurs in the frequency region affected 
by the increase of | Vfi |~^. This shoulder seems to be a 
characteristic of two-spin-waver«rocesses because it has 
been observed in other systemsEEl. 

As the temperature increases, the width of the spec- 
tral function S^^{q, lu) decreases but its height increases. 
The spin wave calculation seems to agree well with the 
MC-SD data for large wavevectors, even for T = 0.6J. 
For small wavevectors and higher temperature, however, 
S'^^(q, oj) from the SCHA calculation is smaller than the 
MC-SD data (Fig. |7|), suggesting that at high tempera- 
tures other processes could be contributing to the dynam- 
ical properties of this system. For systems with easy-axis 
anisotropy one can expect the formation of domains, as 
in the bisip-dimensional Ising model, and, also, localized 
solitonsciJ. In particular, it is usually expectedEHl that 
localized solitons would contribute to the dynamical cor- 
relation function in the a; — > (central peak) region and. 



mostly, for small wavevectors. 



V. CONCLUSIONS 

We have applied a self-consistent harmonic approxi- 
mation to the easy-axis model, obtaining the spinwave 
energies, critical temperature and dynamic correlation 
functions. We also demonstrated how it is possible to 
apply the Wolff cluster Monte Carlo scheme to this easy- 
axis model, by having it act on only the 5^ spin com- 
ponents. For the critical temperature, the SCHA and 
MC results agree favorably over a wide range of easy- 
axis anisotropy, both giving Tc increasing linearly with 
K for K ^ J. The spin-dynamics calculation of dy- 
namic correlation functions shows interesting multiple- 
peak features in ^^^(q, w), that are most easily seen in 
small lattices. These finite-size dynamical features are 
correctly described by the SCHA, especially for T far be- 
low Tc- Similar features should appear in models with 
other symmetries: there are strong evidences that these 
effects were also observed in other simulations of two di- 
mensional easy-plane models.ll3 

All the dynamical calculations discussed in this work 
were performed for anisotropy parameter K = 0.05, 
which corresponds to a transition temperature Tc — 
0.75J. For this anisotropy, two temperatures were an- 
alyzed: T = 0.3J < Tc, and T = 0.6J. We could 
not expect that the spinwave calculation performed here, 
which neglects higher order terms in the spin interactions, 
would reproduce exactly the simulation data. However, 
the agreement for the lowest temperature, T — 0.3J, is 
very good. It is also surprisingly good for T — 0.6J, a 
relatively high temperature, and large wavevectors where 
a lateral peak is seen to develop. At T = 0.6J , for 
small wavevectors and small frequencies, the SCHA func- 
tion for S^^ shows a central peak with height smaller 
than the one obtained from MC-SD simulation. On the 
other hand, the SCHA prediction for the integrated in- 
tensity /^^ for the in-plane correlations lies above the 
MC-SD data for T = 0.6J. These features may suggest 
that other excitations, like localized solitons and domain 
walls, may contribute to the dynamical correlation func- 
tion as the temperature approaches the critical temper- 
ature. It was shown! that the density of these local- 
ized solitons increases exponentially with T as T —^ Tc 
and, then, one should expect that their contribution to 
the dynamics of the system becomes more important for 
temperatures T r^ Tc- To stress this conclusion, we re- 
mark that iS'^^(q, w) obtained by MC-SD simulation for 
T > Tc (not shown here because SCHA cannot be com- 
pared in this temperature regime) does show a central 
peak {uj ^ 0) that increases with T. Its properties will 
be analyzed in a future work. 

We conclude by saying that the two-spin wave cal- 
culation can explain the main features obtained from 
Monte-Carlo-spin-dynamics simulation at very low tem- 



peratures. As T ^ Tc, the comparison between SCHA 
spin-wave calculation and the numerical simulation data 
suggests that other excitations may contribute to the dy- 
namic properties of the model. However, a better under- 
standing concerning the contributions these excitations 
might give to the dynamic spectral functions requires 
some theory which takes into account the existence of 
such objects. To our knowledge, such theory for easy- 
axis anisotropy two-dimensional systems is not available 
in the literature. 
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FIG. 1. The curves correspond to the spin wave energy 
(from (|)) for T = 0.3J (continuous), and T = 0.6J (dashed) 
for K = 0.05J. The circles and stars correspond to the val- 
ues extracted from our numerical simulations; error bars are 
smaller than the symbols. 




FIG. 3. Critical temperature as a function of the 

anisotropy parameter K/J. The symbols correspond to the 
values obtained in our MC calculation (as described in Section 
II B); error bars are smaller than the symbols. 
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FIG. 2. Magnetization as a function of temperature for 
K = 0.05 J. Solid curve is the SCHA theory. Various symbols 
correspond to MC simulation for indicated system sizes; error 
bars are smaller than the symbols. 
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FIG. 4. In-plane integrated intensity I^^ versus wavevec- 
tor, from SCHA (curves) compared with MC-SD (symbols) 
for K = 0.05J, L = 128. 
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FIG. 5. ^^^(q, oj) obtained from numerical simulation for 
' = 0.05J, T 
q= (0.393,0). 
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FIG. 7. S^^{q^,ij) from (continuous line) continuum limit 
(EOl) and from (circles and triangles) numerical simulation 
for K = 0.05J, T = 0.3J, and T = 0.6J, and wavevector 
q= (0.393,0). 
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FIG. 6. ^^^(q,^;) obtained from discrete summation (bs) 
for K = 0.05J, T = 0.3 J, and L = 16, 32, 64, and wavevector 
q= (0.393,0). 
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FIG. 8. ^^^(q,^;) from (continuous line) continuum limit 
( pO[ ) and from (circles and triangles) numerical simulation 
for K = 0.05J, T = 0.3J, and T = 0.6J, and wavevector 
q= (1.473,0). 
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FIG. 9. S^^{(\,ixj) obtained from (continuous line) con- 
tinuum limit (pd) and from (circles and triangles) numerical 
simulation for K = 0.05J, T = 0.3J, and T = 0.6J, and 
wavevector q — (2.503, 0). 



FIG. 11. ^^^(q,^;) obtained from numerical simulation 
(filled circles), discrete summation (empty circles) and from 
continuum limit (line) for K = 0.05J, T = 0.6J, L = 128, 
and wavevector q = (1.03, 0). 
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FIG. 10. S^^{q^,u}) obtained from (continuous line) con- 
tinuum limit (M) and from (circles and triangles) numerical 
simulation for K = 0.05J, T = 0.3J, and T = 0.6J, and 
wavevector q — (1.03, 1.03). 
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FIG. 13. Contours of difference frequency fl{q, k) for 
T = 0.3J and q = (0, 0.393) as function of k. 




FIG. 14. Contours of difference frequency r2(q, k) for 
T = 0.3J and q = (1.03, 1.03) as function of k. 
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